Quantum Phases of Dipolar Bosons in Bilayer Geometry 
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We investigate the quantum phases of hard-core dipolar bosons confined to a square lattice in a 
bilayer geometry. Using exact theoretical techniques, we discuss the many-body effects resulting from 
pairing of particles across layers at finite density, including a novel pair supersolid phase, superfluid 
and solid phases. These results are of direct relevance to experiments with polar molecules and 
atoms with large magnetic dipole moments trapped in optical lattices. 



Recent experimental breakthroughs in the realization 
of ultracold gases of high-spin aligned atoms with large 
dipole moments [IJH], highly excited Rydberg atoms 
[5] , and of ground-state polar molecules [3 H] hold con- 
siderable promise for investigations of many-body quan- 
tum systems where dipolar interactions can become dom- 
inant [9 12 . The anisotropy of dipolar interactions com- 
bined with the possibility to confine particles in low 
dimensional geometries using optical lattices allow for 
study of novel pairing mechanisms and the associated 
quantum phases in a setup where collisional losses are 
suppressed. This is particularly intriguing for the case 
of magnetic atoms, where confinement to lattices with 
spacings as small as 200nm is possible [T3j, which favors 
inter-site dipolar interactions and pairing. 

Pairing of two spin-polarized fermionic dipoles across 
coupled two-dimensional (2D) layers [14] or one- 
dimensional (ID) wires [15] in an optical lattice has 
already lead to the prediction of 2D inter-layer super- 
fluidity [16-21 , analogous to bi-exciton condensation, 
and the ID quantum roughening transition I22j in the 
case of equal number of particles in each layers. Ad- 
ditional exotic phenomena occur for unequal popula- 
tions |23j . where (spin-rotational) symmetry breaking 
can induce, e.g., stable liquids and crystals of compos- 
ite multimers |24j . For bosonic gases in the strongly in- 
teracting regime [25) . emergent parafermionic behavior 
has been demonstrated [261 l2Z] in coupled ID wires. In 
two dimensions, a recent mean-field study in bilayer ge- 
ometry |28j has predicted novel quantum phenomena for 
a model of dipolar bosons on a lattice, including a so- 
called pair supersolid phase. Different from supersolids 
on a single lattice [29TI3T] . the latter implies diagonal 
solid order coexisting with an off-diagonal superfluid or- 
der, both derived from composite pairs of dipoles. The 
experimental observation of this quantum phase and the 
associated pair superfluids and solids would constitute a 
breakthrough for condensed matter in the cold atomic 
and molecular context. Thus, the challenge is now to de- 
termine whether these quantum phases can be realized 
for realistic Hamiltonian representing the microscopic dy- 



namics of strongly interacting dipolar bosons as realized 
in experiments. 

In the present work we study a system consisting of 
hardcore dipolar bosons confined to two neighboring two- 
dimensional (2D) layers of a ID optical lattice. The 
dipole moment of each particle is polarized perpendic- 
ular to the layers, which results in repulsive in-plane 
dipole-dipole interactions. This ensures collisional stabil- 
ity against short-range inelastic collisions in the strongly 
interacting gas. Out-of-plane dipolar interactions are 
dominantly attractive, which favors inter-layer pairing. 
Using exact theoretical techniques based on quantum 
Monte-Carlo methods [32], below we demonstrate that 
this anisotropy and the long-range nature of interactions 
can induce crystallization of the dipolar cloud into a 
charge-density wave for a wide range of trapping parame- 
ters and interactions. Exotic quantum phases such as the 
pair-supersolid (PSS) phase and a pair-superfluid (PSF) 
are achieved under experimentally realistic trapping con- 
ditions. These phases can survive up to temperatures of 
the order of a few nK for a gas of polar molecules or 
strongly magnetic atoms. 

The system we have in mind is described by the single- 
band tight-binding Hamiltonian 

H = — J ^ 4a a ja - ^ V i»;jf} n ^n jyfS 

~y^q n ija . (1) 

Here a,/3 = 1,2 and i,j label the layers and the lattice 
sites in each layer, respectively, while Oj a (aj ) are the 

bosonic creation (annihilation) operators, with aj 2 a = 

and n i a = a\ a a i a . The brackets <> denote summation 
over nearest neighbors only. The first term in Eq. de- 
scribes the kinetic energy with in-plane hopping rate J. 
The second term is the dipole-dipole interaction given 
by Via-jp = C dd (l ~ 3cos 2 6»)/(47r|i Q - j p \ 3 ), where 6 
is the angle between particles at positions i a and jp 
and Cdd = d 2 /e a (Cdd = Mo^ 2 ) f° r electric (magnetic) 
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dipoles of strength d. We denote the repulsive (attrac- 
tive) nearest neighbor intra-layer (inter-layer) interaction 
by V dd = Cdd/(47ra 3 ) (V^ = 2C dd /4wd 3 z ), with a the in- 
plane lattice constant. The interlayer distance is d z . The 
relative strength Vdd/KM can be tuned over a wide range 
of values by changing d z /a. The quantity [i a is the chem- 
ical potential which sets the number of particles in each 
layer. Here we fix /ii = /12, i-c. N\ = N2 

Hamiltonian ([I]) provides a microscopic description for 
the dynamics of, e.g., a gas of RbCs molecules (d « 
1.25 Debye) at low-density n, such that the initial sys- 
tem has no doubly occupied sites [33] . Collisional stabil- 
ity is ensured for rT 1 ! 2 ^> (d 2 /huj^) 1 / 3 ~ 130nm with 
cjj_ ~ 100kHz the frequency of transverse confinement 
provided by the in-plane optical lattice [3I| ■ In addition, 
the choice d 2 /d\ < Vo avoids interaction- induced inter- 
layer tunneling, with Vo the depth of the optical potential 
in the transverse direction. Model ([!]) can also be used 
to describe the dynamics of a gas of strongly magnetic 
dipolar atoms, such as Dy (d — 10/is). In this case the 
conservative estimate above for collisional stability is sat- 
isfied for lu± ~ 1 kHz. 

In the following, we present exact theoretical results 
based on path integral Quantum Monte Carlo simula- 
tions using a two-worm algorithm |35j which allows for 
efficient sampling of paired phases. We have performed 
simulations of L x L = N sites square lattices with 
L — 8, 12, 16, 20 and 24. For computational convenience, 
we have set the dipole-dipole interaction cutoff to the 
third nearest neighbor and have checked that using 
a larger cutoff did not change the simulation results 
within errorbars. Lower cutoff values do not allow for 
stabilization of, e.g., supersolid phases, see below. In 
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FIG. 1. (Color online) Phase diagram of Hamiltonian |T]) 
as a function of Vdd /J and particle density n, computed via 
QMC simulations, for an interlayer distance d z /a — 0.36 (see 
text). CB: Checkerboard solid; PSS: Pair supersolid; PSF: 
Paired superfluid; 2SF: independent superfluids. The phase 
boundaries in the dashed region are not resolved. 
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TABLE I. Quantum phases of Fig. [T] and corresponding or- 
der parameters: structure factor S(n,ir); single-particle con- 
densate ipf = (a i>a ) in each layer a; pair-condensate order 
parameter ^ = {a ita a it p), with a / /3. (See text) 



the following we choose d z /a — 0.36 and d z ~ 200 nm, 
which is experimentally feasible with, e.g., Cr or Dy 
atoms [H [3J . We show below that this choice allows one 
to access a parameter regime where particles on different 
layers can pair up to form a composite object. Below, 
we first discuss the phase diagram, and then discuss in 
more details the various phases. 

The phase diagram of Eq. ([I]) at temperature T = 
is shown in Fig. [I] as a function of Vdd/ J and the 
density n, in the parameter regime 0.31 > Vdd/ J > 0.2 
and 0.1 < n < 0.9, with d z /a = 0.36. We expect 
this phase diagram to be representative of situations 
with d z /a <Si 2, where interlayer pairing is favored (see 
figure [2] below). 

At half-filling n = 0.5, an incompressible checkerboard 
solid of pairs (PCB) is stabilized for sufficiently large val- 
ues of Vdd/ J- Similar to the conventional checkerboard 
phase present in single-layers (3D], here atoms in each 
layer occupy every other site of the lattice, due to in- 
plane dipolar repulsion. The checkerboard order is char- 
acterized by a finite value of the static structure factor 
5(k) at the reciprocal lattice vector k = (ir,ir), with 

S( k ) = jf E ex Pt ik ( r " (2) 

r,r/ 

and the system displays zero superfluidity. We find that 
in the PCB phase atoms across the layers are strongly 
paired due to attractive interlayer interactions. As a 
result, the position of the two checkerboard solids is 
strongly correlated, i.e., they sit on top of each other. 
The system can be thus envisioned as a solid of pairs 
[3D], with an effective mass m e s ~ J 2 /(2V d J d + zV d d), 
where z is the coordination number. The PCB solid is 
stabilized at (much) lower values of Vdd/ J compared to 
the case of checkerboard solids in a single layer [3D], in 
analogy with what found in [36j . This is due to the higher 
effective mass of the pairs. A similar robustness of this 
phase is also found for melting at finite temperature. 

Upon doping the PCB solid with extra particles or 
holes, a so-called pair-supersolid (PSS) phase is immedi- 
ately stabilized. The latter displays both diagonal long 
range order with S^tt, 7t) 7^ 0, off-diagonal long-range or- 
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der associated with a non-vanishing value of the pair- 
condensate order parameter = (a% >a ai p) 7^ (with 
a 7^ 0), and an associated finite superfluid stiffness for 
pairs (see below). The single-particle condensate order 
parameter ipf — (a,i >a ) — is instead zero. The exis- 
tence of off-diagonal order is consistent with a picture of 
delocalized defects 29, 37 , which here correspond to cor- 
related pairs of holes or extra particles across the layers. 
The PSS phase forms a lobe structure in the (Vdd/J—n)- 
plane, around the PCB line. Away from the tip of the 
lobe, we find that by varying n at constant Vdd/ J the 
PSS loses its diagonal long range order by melting into a 
pair superfluid phase (PSF), via an Ising-type transition 
(red continuous line). The PSF phase, with ^ 7^ and 
ijjf = 0, is destroyed in favor of a 2SF (a phase with in- 
dependent, though correlated, superfluids on each layer) 
for smaller values of Vdd/ J- In particular, we notice that 
a tiny PSF-region should persist in between the PSS and 
2SF phases even close to the tip of the PSS-lobe, however 
this is within errorbars for Vdd/ J ^5 0.2. Exactly at filling 
n = 0.5, our results are consistent with a direct PCB-2SF 
transition, as discussed below, with no intermediate PSS 
phase. In particular we find no evidence of, e.g., possible 
micro-emulsion phases [3TJ [38] , within errorbars. 

Finally, we notice that a host of other phases are 
present in the general phase diagram for two layers. In 
particular, we find that for stronger values of Vdd/ J > 0.3 
the system displays a sequence of incompressible phases 
at various rational fillings of the lattice, similar to the 
so-called Devil's staircase found in the case of a single 
layer. We also expect novel PSS phases to appear around 
lobes at, e.g., filling n = 0.25, in analogy with Ref. [3"U] . 
In addition, independent solids as well as supersolid 
phases can be achieved by increasing the layer distance, 
while mixtures of solid and superfluid phases can be 
stabilized by modifying the relative particle density in 
the two layers. The discussion of some of these phases 
is however outside of the scope of the present work. In 
the remainder of the paper we discuss in more detail the 
various phases and their transitions at zero and finite 
temperature around n = 0.5. 

Stability of the PCB phase: As discussed above, the 
PCB phase at n — 0.5 is characterized by a finite value 
of the order parameter S(tt,t:) and no off-diagonal or- 
der. The latter is associated with superfluidity in a (2+1) 
dimensional interacting system, which can be measured 
straightforwardly within Monte-Carlo (see below) . In ad- 
dition, within the PCB phase inter-layer dipolar attrac- 
tion strongly correlates the positions of particles in the 
two layers. 

The stability of the PCB phase with respect to intra- 
plane interactions as well as inter-layer distance d z /a 
at zero temperature is analyzed in Fig. [2j There, we 
numerically determine the minimum dipolar interaction 
strength Vdd/ J required to stabilize the PCB phase at a 



given d z Ja. In order to establish whether the solid phase 
is paired we have performed several simulations with dif- 
ferent initial conditions for each set of parameters and 
observed whether the equilibrium configuration was de- 
pendent on the initial choice or not. The figure shows 
that a PCB phase is stabilized for d z /a < 2 and suffi- 
ciently large Vdd/ J (continuous line). In this parameter 
regime, the system above (below) the continuous line is 
a PCB (2SF) phase, respectively, that is, the continuous 
line visualizes the shift of the PCB-2SF transition point 
of Fig.[l]as a function of d z /a. Instead, for d z /a > 2 and 
large enough interactions the insulating phase above the 
(dotted) line corresponds to two independent checker- 
board phases (2CB). This points to the possible pres- 
ence of a tri-critical point in the phase diagram around 
d z /a « 2. We have confirmed that the computed transi- 
tion points are independent of the interaction cutoff that 
we use, within our errorbars, and should be thus quanti- 
tatively relevant to experiments. 

In the following we focus on d z /a = 0.36 to satisfy 
Vtid ~ 10 J in the vicinity of the tip of the lobe. We find 
that this choice ensures pairing at n ~ 0.5 (in the vicin- 
ity of the PCB phase) while keeping Vdd relatively low. 
This corresponds to experimentally optimal conditions to 
observe PSS phase: a lower effective mass of pairs m e g 
results in a larger superfluid density which in turn results 
in higher critical temperatures (see also below). 
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FIG. 2. (Color online) The plot of minimum Vdd/ J needed to 
stabilize the CB phase as a function of d z /a. Once the layers 
are separated by d z /a > 2 they behave as independent layers. 

Pair supersolid phase: Figure [T] shows that a PSS lobe 
is immediately formed by doping the PCB solid with ei- 
ther vacancies (holes) or interstitials (extra particles). 
The hard-core constraint of Eq. ([!]) ensures particle-hole 
symmetry, and thus reflection symmetry of the lobe, 
around n = 0.5. 

We characterize this pair supersolid phase in Fig. [3j for 
a specific choice of interaction strength Vdd/ J = 0.238. In 
the figure, the order parameter for the diagonal checker- 
board solid order S(ir, tt) (continuous lines) and the su- 
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perfluid stiffness of pairs ppss are plotted as a function 
of n. The quantity p PSS = T(W 2 )/dL d - 2 [39] is directly 
related to a pair condensate, and can be calculated within 
quantum Monte-Carlo, with W = W\ + W% the sum of 
winding numbers in layer 1 and 2. The figure shows 
that for an extended range of densities, both the static 
structure factor and the pair superfluid stiffness are fi- 
nite and system size independent, showing the existence 
of a stable supersolid phase in the lobe region. We note 
that, due to pairing across the layers, in the PSS phase 
the fluctuation of difference in winding numbers is zero 
((W 1 -W 2 ) 2 ). 

Superfluid phases: As the system is doped further, 
the PSS disappears in favor of a PSF phase. The 
latter displays pair-induced off-diagonal long range 
order, only [see, e.g., Table 1]. We find that the 
PSS-PSF transition is of the Ising type universality 
class in (2+l)-dimensions, analogous to the case of a 
single-layer |30) . Critical points are determined using 
finite size scaling for the static structure factor with 
scaling coefficients 2j3/v = 1.0366 gD] (see Fig. ||b) 
for the specific choice Vdd = 0.238 J). In the figure the 
scaled quantity S(ir, tt)L 1 0366 is plotted as a function of 
n, and the crossing of the curves at n CI — 0.573 ± 0.002 
corresponds to the quantum critical point where the 
finite size effects disappear [see also panel (a)]. 

We find in general that by lowering the interaction 
strength Vdd/ J at constant n from the PSF phase, the 
system finally develops into two independent superflu- 
ids (2SF) with a finite value of the single-component 
condensate order parameters, tpf = ^ 0, via a sec- 
ond order phase transition in the (2+1) XY universality 
class. The transition points between the PSF and 2SF 
phases in Fig. [I] are calculated using finite size scaling of 
((Wi - W 2 ) 2 ). The latter quantity is zero inside the PSF 
phase in the thermodynamic limit due to pairing across 
the layers, while it has a finite value in the 2SF phase. 
We note that the pair order parameter in the 2SF phase 
is instead trivially non-zero, \& ^ (see also Table 1). 

The phase diagram in Fig. [I] shows that the boundary 
of the PSF-2SF transition shifts downward approxi- 
mately linearly in the (Vdd/ 'J ~ fi)-plane, as the density 
becomes sufficiently smaller or larger than n — 0.5. 
This is easily understood in the limit of very small 
densities, by noting that inter-plane dipole-dipole 
interactions always favor the existence of a two-body 
bound state, even for an arbitrarily small interaction 
strength. However, we find that many-body effects 
result in a threshold for the formation of pairs at 
finite density, where the magnitude of the interaction 
strength required to stabilize pairing increases with n. 
This is explained by noting that, in the limit of low 
density, d z y/n <C 1, PSF phase is composed of weakly 
interacting superfluid dimers. As the density is increased 
exchanges between dimers are favored. This destabilizes 




FIG. 3. (Color online) (a) Structure factor S(n, n) (solid 
lines, left y-axis) and superfluid stiffness p s (dashed lines, 
right y-axis) in the PSS phase for L — 8 (black squares), 
12 (red circles), 16 (blue triangles) and 20 (green diamonds) 
at T/J = 1/(1. 5L) shown using black squares, red circles, 
blue triangles and green diamonds respectively. The PSS- 
PSF transition point is at Vdd = 0.238 J. (b) Scaled structure 



factor S(jv, ty)L 
and 20. 



n with 2/3/i/ = 1.0366 for L = 8, 12, 16 



the dimers, inducing the transition to two independent 
superfluids in the 2SF phase. Eventually, the presence 
of diagonal order near n = 0.5 forces the PSF-2SF line 
to bend down, deviating from the linear dependence on n. 

We gain further insight into the structure of correla- 
tions in the condensed phases by studying the following 
four-point correlation function: 



fji = (^i,iip2,iWi 



i4,l) 



(3) 



Here i,j, I refer to sites, 1, 2 refer to layers, and () denotes 
a quantum and thermal average as well as site averaging 
over i. In the presence of pair superfluidity, one expects 
this correlation function to be short ranged with respect 
to rji = \rj — r{\, and simultaneously long ranged with 



respect to ru 



n\ and n 



In the 2SF 



phase, instead, fji is obviously long ranged with respect 
to ru and , but it is independent of Tji . 

Figure [4] shows fji (normalized to unity: 
/o°° fjld r jl = 1) as a function of rji for the PSS 
(green triangles, n = 0.48, V dd = 0.25 J), PSF (red dots, 
n = 0.40, Vdd = 0.25J) and 2SF phases (blue squares, 
n = 0.30, Vdd = 0.18J). As expected, we find that fji 
is independent of in the 2SF phase, where pairing is 
absent, while it is peaked at ru = both in the PSS 
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FIG. 4. (Color online) Four-point correlation function fji of 
Eq. f3) as a function of Tji for the 2SF (blue squares, n = 
0.30, V Ad = 0.18 J), PSF (red dots, n = 0.40, V dd = 0.25 J) 
and PSS (green triangles, n = 0.48, Vdd = 0.25J) phases. 
The dashed (dotted) line is the exponential fit, foe~ Tjl ^ a , to 
the PSF (PSS) histogram, where £o can be interpreted as the 
extent of the pair wavefunction (see text). The inset shows 
£o across the PSF-PSS phase boundary. 



and PSF phase. The figure shows that an exponential 
ansatz of the form foe~ Tjl ^o fits quite well the large-r^ 
behavior of fji in these latter phases, and is essentially 
exact for all rji in the PSS phase with £o — 1.63a. 
Here £o can be interpreted as the spread of the pair 
wavefunction, and is obtained from Fig. [4] by fitting the 
tail of fji , as obtained numerically. The inset in Fig. [4] 
shows £ as a function of n, as the PSF-PSS phase 
boundary is crossed. The pair wavefunction is shown to 
be considerably more tightly bound in the PSS phase 
than in the PSF phase. The abrupt drop in £o locates 
precisely the transition point. 



Finite temperature: We have studied the robustness 
of the quantum phases described above against thermal 
fluctuations. As expected for two-dimensional systems, 
we find in general that superfluidity in the PSS, PSF 
and 2SF phases disappears at finite temperature T via a 
Kosterlitz-Thouless (KT) type [UJ transition. Diagonal 
long range order in the PCB and PSS phases is instead 
lost via a two-dimensional Ising-type transition. We have 
found that, when present, pairing still exists at the tran- 
sition points, suggesting that the temperatures required 
for breaking pairs are higher than the critical tempera- 
tures measured here. 

Figure [5] shows one example for the SF- normal transi- 
tion in the 2SF phase. We plot p s vs. T/J at Vdd/ J = 
0.20 and n = 0.3 for different system sizes. The inset 



shows the finite size scaling procedure [12] used to de- 
termine the critical temperature. We find Tkt,2SF = 
Trh 2 Ps (T KT ) _ Q 255J. For the PSF-normal transition we 
find T kt ,psf ~ 0.08J at n = 0.3 and V dd /J = 0.25. 
The lower KT transition temperature compared to the 
2SF-normal transition is due to a larger effective mass of 
the pairs, i.e. lower effective hopping, which results in 
a suppression of particle delocalization and consequently 
smaller p s . The disappearance of the PSS phase pro- 
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FIG. 5. (Color online) Superfluid stiffness p 3 as a function 
of temperature T/J, at Vdd = 0.20J and n — 0.30, corre- 
sponding to 2SF phase at L — 12, 16,20 and 24 shown us- 
ing black squares, red circles, blue triangles and green di- 
amonds respectively. As temperature is increased the 2SF 
phase undergoes KT phase transition at critical temperature 
Tkt,2Sf ~ 0.255 J, indicated by an arrow. The inset shows 
finite size scaling [?2] where the dashed line is a linear fit of 
our simulation results (points). 



ceeds in two successive stages. At Tkt,pss the PSS phase 
melts into a liquid-like phase reminiscent of a liquid crys- 
tal, with p s = and S(jr, tt) ^ 0. Upon further increasing 
the temperature S(ir, n) becomes zero at a critical tem- 
perature T c through an Ising-type transition {2f3 /v = 1/4 
in 2D). For example, we find Txt.pss ~ 0.06 J and 
T c w 0.3J for n = 0.48, V dd = 0.25J. Similar T c val- 
ues are found for the critical temperature of the melting 
of the PCB phase into a featureless normal fluid, e.g., 
T c sa 0.35J for V d d = 0.25J. Clearly, for larger interac- 
tion strengths, i.e. away from the tip of the lobe, transi- 
tion temperatures will increase. 

Experimental estimates: Based on our results we es- 
timate under which experimental conditions the phases 
described can be observed. For example, with a gas of Dy 
(d = lOps) a choice of lattice parameters a — 500 nm, 
d z = 200 nm, J = 50/iHz results in V dd /J - 0.21 
which stabilizes the PCB phase. In the case of Er2 Fes- 
hbach molecules [?31lll](d = Up B ) with a = 400 nm, 
d z = 200 nm, J = 100 kHz the PCB phase is stabilized 
at Vdd/J ~ 0.4. In both cases the PCB phase can be 
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observed at nk temperatures. 

Using RbCs (d = 0.3D) and typical trapping parame- 
ters a = 500 nm, d z = 300 nm and J = 150 hUz we find 
Vddl J = 0.7, which is large enough to stabilize the PCB. 
The latter survives up to T^ CB ~ 4 nK. By doping away 
from filling factor n = 0.5 the PSS phase can be reached 
with a KT transition temperature for PSF-normal tran- 
sition of the order of nK. 

In conclusion, we have studied the quantum phases of 
dipolar bosons in a bilayer lattice geometry described 
by the microsopic Hamiltonian Eq. ([I]) for hard-core 
particles, in a situation where the number of particles 
in each layer is the same. Relevant to experiments with 
polar molecules and magnetic atoms, we have estab- 
lished under which conditions pairing for two particles 
is stabilized across the layers. Our zero temperature 
study indicates that the system displays a rich ground 
state phase diagram including a novel pair-supersolid 
phase for hard-core dipolar bosons, in addition to pair 
superfluid and checkerboard-like solid phases. Our 
finite temperature results indicate that these phases are 
experimentally observable at temperatures of the order 
of nK. A four-body correlation function connected with 
the spread of the pair wave-function can be used to 
characterize these phases and their transitions. Future 
work will include the extension of similar quantum 
Monte Carlo studies to multilayer geometries as well as 
to systems with population imbalance in the layers. 
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